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1 Abstract 

We perform numerical simulations of the 2-d antiferromagnetic quantum 
Heisenberg model using an efficient cluster algorithm. Comparing the finite 
size and finite temperature effects of various quantities with recent results from 
$h ' chiral perturbation theory we are able to determine the low energy parameters 

of the system very precisely. We find eo = —0.6693(1) J /a 2 for the ground 
state energy density, M s = 0.3074(4)/a 2 for the staggered magnetization, 
he = 1.68(l)Ja for the spin wave velocity and p s = 0.186(4)J for the spin 
stiffness. Our results agree with experimental data for the undoped precursor 
insulators of high-T c superconductors. 
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The first high-T c superconductor to be discovered was La2- x Ba x .Cu04 with x ~ 
0.15 which has a layered structure with 2-d copper-oxygen planes. The copper 
ions are located at the sites of a quadratic lattice with lattice spacing a = 3.79 A. 
The undoped material La2Cu04 is an insulator, however, with strong antiferro mag- 
netic interactions within the copper-oxygen planes between electron spins localized 
at the copper ions. The couplings between different layers are extremely weak. 
Experimentally one observes long range antiferromagnetic order, i.e. a spontaneous 
staggered magnetization M. s arises, which breaks the 0(3) spin rotational symmetry 
down to 0(2). The low energy excitations of the system are spinwaves (the so-called 
magnons) which are the Goldstone bosons of the spontaneously broken 0(3) sym- 
metry. The physical situation can be modeled by the 2-d Heisenberg quantum spin 
system with Hamiltonian 

H — J 2_, S x ■ S x+ p,, (1) 

where S x = ^a x is a spin ^ operator {a x are Pauli matrices) located at the point x 
of a 2-d quadratic lattice with lattice spacing a. The interaction is between nearest 
neighbors {fx is the unit vector in //-direction) and J > is the antiferromagnetic 
exchange coupling. The question arises how well this model describes the physics of 
the copper-oxygen planes in La2Cu04, in particular how it compares quantitatively 
with experimental results. 

Here we concentrate on the calculation of the low energy parameters of the 
model, which determine the dynamics of the Goldstone bosons. These are the 
staggered magnetization Ai s , the spinwave velocity he and the spin stiffness p s . 
Based on symmetry considerations chiral perturbation theory makes very strong 
predictions for the magnon dynamics, containing the low energy parameters as the 
only unknown constants. Recently, Hasenfratz and Niedermayer have worked out 
the chiral perturbation theory for the antiferromagnet in great detail up to two- 
loop order 0. Lower order results had been obtained before by Fisher || and by 
Neuberger and Ziman || . Here we only quote the results of ref . that are essential 
for our study. We consider the system at finite temperature T and in a finite spatial 
volume of size L x L with periodic boundary conditions, with I 3 = hc/TL such that 
I is of order 1. For small enough temperatures T <C 27rp s and large enough volumes 
hc/L <C 2iip s Hasenfratz and Niedermayer obtained the following results: for the 
internal energy density 



e{T,L) = e -^-{l + 4(3 {l) '"' 



+ ■■■}, (2) 



3L 2 { dr " ' p s Ll 
where e is the ground state energy density; for the staggered susceptibility 

Xs{T, L) = ^ J 1 + + \mf + 3&(l) +...}; (3) 
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and for the uniform susceptibility 



2p s J l he ~ l f he 



2 r 



x (T,L) = - 7 £z i+o-tt /MO + o 



3(/ic) 2 I 3 p s Ll 3 \p s Ll 



, (4) 

The functions /3;(7) and ^(Z) are shape coefficients which depend only on / and 
which are described in detail in ref. @ . We will use these results of chiral perturbation 
theory to determine the unknown low energy parameters eo, M. S) he and p s from a 
fit of e(T, L), Xs{T, L) and x(T, L) to numerical Monte Carlo data. This method has 
been used before for classical spin models and for relativistic quantum field theories 



First we decompose the Hamiltonian into H = Hi + H 2 + H 3 + H^ with 
Hi = J ^2 &x- S x+ i, H 2 = J ^2 x S x+ 2, 

x=(2m,n) x=(m,2n) 

H 3 = J ^2 Sx 4 S x +b H 4 = J ^2 Sic' S x+ 2, (5) 

i=(2m+l,n) x=(m,2n+l) 

and we use the Suzuki- Trotter formula for the partition function 

Z = Trexp(-(3H) = lim Tr[exp(-e/3ifi) exp(-e(3H 2 ) exp(-e(3H 3 ) exp(-e(3H 4 )} N , 

N— >oo 

(6) 

where (3 = 1/T is the inverse temperature and e = 1/N determines the lattice 
spacing in the euclidean time direction. By inserting complete sets of eigenstates 
|1) and |-1) of a x between the factors exp(— e(3H,j) we map the 2-d quantum spin 
system to a 3-d induced classical system of Ising-like variables s(x, t) = ±1 (t labels 
the euclidean time slice) 

Z = U E eM-s) (7) 

x,t s(x,t)=±l 

with an action 

S = J2 S[s(x,t),s(x + l,t),s(x,t + l),s(x + l,t+ 1)] 

x=(2m,n) ,t=Ap 

+ S[s(x,t),s(x + 2,t),s(x,t + l),s(x + 2,t + l)} 

x=(m,2n) ,t=4p+l 

+ J2 S[s(x,t),s(x + l,t),s(x,t + l),s(x+l,t+l)) 

x=(2m+l,n),t=4p+2 

+ S[s(x,t),s(x + 2,t),s(x,t + l),s(x + 2,t + l)}. (8) 

x=(m,2n+l),t=4p+3 

The classical spins interact with each other via four-spin couplings S[s(x, t), s(x + 
p, t), s(x, t + 1), s(x + p,t + 1)] associated with time-like plaquettes. Up to a trivial 
additive constant one has £[1,1,1,1] = S[— 1, — 1, — 1, — 1] = 0, S[l, — 1, 1, — 1] = 
£[-1,1,-1,1] = -log[|(exp(e/3J) + 1)] and £[1,-1,-1,1] = £[-1,1,1,-1] = 



— log[|(exp(e/3J) — 1)]. All other action values are infinite. This causes problems 
in numerical simulations because many spin configurations are forbidden and the 
updating must respect several constraints. In a previous paper we have introduced 
blockspins |J to resolve the constraints. For the 1-d antiferromagnetic spin chain 
the blockspin model is not frustrated and the use of a blockspin cluster algorithm 
eliminates critical slowing down. In two dimensions, however, frustration causes 
severe problems. Recently, Evertz, Lana and Marcu have developed loop cluster 
algorithms for vertex models, which can also be applied to quantum spin systems. 
The algorithm constructs closed loops of spins and flips them simultaneously. The 
loop cluster algorithm does not suffer from frustration but it may suffer from so- 
called freezing. Freezing occurs when a loop branches out many times and fills a 
large fraction of the whole volume. We find that freezing does not arise for the 
Heisenberg antiferromagnet. This is essential for the success of our numerical study. 

The algorithm constructs loops by first selecting a starting point (x, t) at random. 
The spin s(x,t) participates in two plaquette interactions, one at euclidean times 
before and one at euclidean times after t. When s(x, t) = 1 we consider the plaquette 
interaction at the later time, and for s(x,t) = —1 we consider the interaction at the 
earlier time. The corresponding plaquette configuration is characterized by the spin 
orientations at the four corners. One of the corners will be the next point on the loop. 
For configurations C\ = [1, 1, 1, 1] or [—1, — 1, —1, —1] the next point is the time-like 
nearest neighbor of (x,t) on the plaquette. For configurations C 2 = [1, — 1, 1, —1] or 
[— 1, 1, —1, 1] the next point on the loop is with probability p = 2/(exp(e/3J) + 1) 
the time-like nearest neighbor, and with probability 1 — p the space-like nearest 
neighbor of (x,t). Finally, for configurations C3 = [1, —1, —1, 1] or [—1, 1, 1, —1] the 
next point on the loop is the space-like nearest neighbor of (x,t). Once the next 
point on the loop is determined the process is repeated until the loop closes. Then all 
spins on the loop are flipped simultaneously. The algorithm obeys detailed balance, 
p{Ci)w{Ci -> Cj) = p{C j )w{C j -> d), where p(Ci) = 1, p(C 2 ) = |(exp(e / SJ) + 1), 
p(Cs) = |(exp(e/3J) — 1) and w(Ci — ► Cj) is the transition probability to go from a 
plaquette configuration Cj to Cj. Indeed one has 

p(CiMCi - C 2 ) = 1 = -p = p(C 2 )w{C 2 - d), 

p 

p(C 2 )w(C 2 - C 3 ) = -(l-P) = i(e*p(e/3J) - 1) = p(C 3 )w(C 3 - C 2 ). (9) 
p 2 

In our construction a loop cannot branch out and hence freezing does not arise. 
Cluster algorithms offer the possibility to use improved estimators which reduce the 
variance of different observables. For example, the uniform susceptibility can be ex- 
pressed as x a2 J — im^c/ |C|)> where AN is the number of points in the euclidean 
time direction, \C\ = J2( x ,t)ec 1 i s ^ ne s i ze °f the loop C and M c = \ H( x ,t)ec s ( x > t) * s 
the loop magnetization. It is interesting to note that clusters with nonzero magneti- 
zation must wrap around the lattice in the euclidean time direction. Small clusters 
which do not wrap around the lattice have Mc = 0. Similarly, one can define im- 
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Table 1: Numerical data for e, Xs and \. 
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256 


-0.678(1) 


9.67(3) 
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256 


-0.673(1) 


16.08(5) 


0.0514(3) 
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10 


256 


-0.672(1) 


23.73(7) 


0.0527(3) 
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-0.672(1) 
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0.0528(3) 
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256 
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0.0535(3) 
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10 


12 
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-0.673(1) 
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0.0460(3) 
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-0.670(1) 
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-0.673(1) 


103.1(4) 
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0.0111(3) 


15 
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-0.675(1) 


34.8(1) 


0.0287(3) 


15 
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-0.674(1) 


57.9(2) 


0.0385(3) 
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-0.674(1) 
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-0.672(1) 
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768 
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0.0451(3) 


15 
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-0.670(1) 
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0.0457(3) 


15 


20 


768 


-0.671(1) 


227.6(8) 


0.0457(3) 



proved estimators for the staggered susceptibility Xs and for the internal energy 
density e. 

Some results of our numerical simulations are collected in table 1. We have per- 
formed measurements for three inverse temperatures f3J = 5, 10, 15 and for different 
spatial sizes L/a = 6,8, ...,20. We have always performed 10000 loop updates for 
equilibration followed by 100000 measurements using the improved estimators. The 
autocorrelation times of the loop cluster algorithm are at most a few sweeps, and 
we see no indication of critical slowing down. With standard local algorithms it 
would be impossible to reach temperatures as low as the ones we use here, because 
of severe problems with slowing down. In table 1 the lattice spacing has been fixed 
to e(3J = ^. We have also performed runs on coarser lattices with e(3J = ^ and 
J^. This allows us to extrapolate our data to the euclidean time continuum limit 
e — > 0. After the extrapolation we fit the results to the above expressions from chiral 
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perturbation theory. The data for e, Xs an d X are an fitted simultaneously. Our 
best fit with x 2 /dof = 1.4 is shown in fig 1. The finite size and finite temperature 
effects of the internal energy density depicted in fig. la are very small (of the order 
of our statistical errors), while the effects on the susceptibilities are much larger. 
For low temperature and small volume some data have been excluded from the fit 
because for them I is not of order 1. The fit gives the following values for the low 
energy parameters 

e = -0.6693(1) J/a 2 , M s = 0.3074(4)/a 2 , he = 1.68(1) Ja, p s = 0.186(4) J. (10) 

To our knowledge this is the most accurate determination of these zero temper- 
ature and infinite volume properties from a simulation of the partition function 
at finite temperature and finite volume. The result for the ground state energy 
density agrees with different zero temperature Monte Carlo calculations || which 
yield e = -0.6692(1) J/a 2 . Our results are consistent with an analytic expansions 
around the Ising limit || which gives e = 0.6693(1) J/a 2 and A4 S = 0.307(l)/a 2 , 
but not consistent with a recent large scale numerical study using a standard lo- 
cal algorithm which obtained p s = 0.199(2)J. Finally, we compare our results 
with experimental data. Using inelastic neutron scattering the spin wave velocity 
he = 0.85(3)eVA has been measured [|ll]], while an analysis of Raman scattering 
data yields J = 0.128(6)eV = 1480(70)K. Using this together with a = 3.79A 



the experiments on La2Cu04 obtain he = 1.75(9)Ja. This is consistent with our 
result for the Heisenberg antiferromagnet. Using the experimental values for the 
spinwave velocity and for the lattice spacing we obtain an independent estimate of 
the exchange coupling in La 2 Cu0 4 

J = 0.133(5)eV = 1540(60)K. (11) 

The agreement between our numerical results and the predictions of chiral pertur- 
bation theory confirms that the Heisenberg model has long range antiferromagnetic 
order, and that its low energy dynamics is dominated by magnons. A precise de- 
termination of the low energy parameters that determine the magnon physics was 
possible only because the loop cluster algorithm is very efficient also at low temper- 
atures. Recently, a loop cluster algorithm has been constructed for lattice fermion 



systems [13]. This raises hopes that numerical investigations of similar accuracy 
become feasible for the Hubbard model, and hence for high-T c superconductors like 
La2- x Ba a; Cu04. 

We are indebted to P. Hasenfratz and F. Niedermayer for making their results 
available to us prior to publication. We also like to thank them for many interesting 
discussions about quantum antiferromagnets. 
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Figure Caption 

Fig.l: The fit of the Monte-Carlo data for the internal energy ea 2 /J (a), the stag- 
gered susceptibility Xs^J/L 2 (b), and the uniform susceptibility x a2 J ( c )- The 
dots, squares and triangles are the Monte-Carlo data for (3J = 5, 10 and 15 re- 
spectively. The corresponding fit functions are represented by the solid, dashed and 
dotted curves. 
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